function result = phi0 (r, r_, theta)
%PHI0 This function gets 3 arguments (r, r', theta) and returns phi0(r, r', theta)
%   This function implements the following formula:
%   phi0(r, r_, theta) = q * integral(0 => 1)[A(sqrt(r'^2(1-alpha)^2 + 
%      alpha^2*r^2 + 2*alpha*(1-alpha)rr'cos(theta)))d-alpha]

% Use global parameter:
global q;
global h_alpha;

result = cachedPhi0(r, r_, theta);

if isnan(result)
    result = q .* trpzInt(@(alpha)int_alpha(alpha, r, r_, theta), 0, 1, h_alpha);  
end

end
